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ABSTRACT 

We present an ultrafast opacity calculator that we name HELIOS-K. It takes a line list as an input, com¬ 
putes the shape of each spectral line and provides an option for grouping an enormous number of lines into a 
manageable number of bins. We implement a combination of Algorithm 916 and Gauss-Hermite quadrature to 
compute the Voigt profile, write the code in CUDA and optimise the computation for graphics processing units 
(GPUs). We restate the theory of the k-distribution method and use it to reduce ^ 10^-10® lines to ^ 10-10'' 
wavenumber bins, which may then be used for radiative transfer, atmospheric retrieval and general circulation 
models. The choice of line-wing cutoff for the Voigt profile is a significant source of error and affects the value 
of the computed flux by ^ 10%. This is an outstanding physical (rather than computational) problem, due to 
our incomplete knowledge of pressure broadening of spectral lines in the far line wings. We emphasize that this 
problem remains regardless of whether one performs line-by-line calculations or uses the k-distribution method 
and affects all calculations of exoplanetary atmospheres requiring the use of wavelength-dependent opacities. 
We elucidate the correlated-k approximation and demonstrate that it applies equally to inhomogeneous atmo¬ 
spheres with a single atomic/molecular species or homogeneous atmospheres with multiple species. Using a 
NVIDIA K20 GPU, HELIOS-K is capable of computing an opacity function with ^ 10® spectral lines in ^ 1 
second and is publicly available as part of the Exoclimes Simulation Platform (ESP; www. exoclime . org). 

Subject headings: radiative transfer — planets and satellites; atmospheres — methods; numerical 


1. INTRODUCTION 

1.1. The million- to billion-line radiative transfer challenge 

Measuring the spectra of exoplanetary atmospheres gives 
us a wind ow into their thermal structure and c hemical com¬ 
positi o ns (iBrowrj 1200 It iBurrows e t al.l 1200 11; ICharbonneaul 
2009t IS eager & Demina 120101 ; iMadhusudhan et al.l 120141 ; 
iHeng & Showmaniboish . A crucial bridge between obser¬ 
vation and inference is the use of theoretical models of atmo¬ 
spheric radiation, both in the form of “forward models” that 
adopt a set of fixed assumptions (e.g., solar composition) and 
retrieval models that attempt to invert for various properties 
from the data. In both families of models, one needs to com¬ 
pute synthetic spectra, which in turn requires the computation 
of the opacity function of the atmosphere. 

To achieve a high degree of accuracy, it is desirable to per¬ 
form “line-by-line” calculations, where every spectral line in 
the range of wavelengths considered, for a given molecule 
(e.g., water), is directly included either in the process of solv¬ 
ing for radiative equilibrium (in forward models) or a multi¬ 
parameter search for an optimal solution based on a compar¬ 
ison to data (in retrieval models). Such an approach may 
be readily adopted at low temperatures, but at the high tem¬ 
peratures (^ 800—3000 K) of the exoplanetary atmospheres 
currently amenable to characterisation by astronomy, it be¬ 
comes infeasible as the number of spectral lines involved in¬ 
creases by orders of magnitude. Eor example, the HITRAN 
database lists ^ 10® lines for the water molecule, but is only 
valid up till temperatures of about 800 K. At higher tempera¬ 
tures, millions of weak lines become important and the total 
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number of lines involved increases to ~ 10®; the HITEMP 
database needs to be used instead. Line-by-line calculations 
become expensive or even prohibitive as one attempts to ex¬ 
plore the broad parameter space occupied by exoplanetary at¬ 
mospheres. Eurthermore, in studies where line-by-line calcu¬ 
lations are claimed, it is not always clear that sufficient res¬ 
olution has been devoted to computing the > 10® lines of 
the opacity function for hot exoplanetary atmospheres. As 
different combinations of molecules, temperature and pres¬ 
sure are considered, the problem becomes computationally 
intractable. 

1.2. The method of k-distribution tables 

In the Earth and planetary sciences, a well-worn strat¬ 
egy for dealing with an enormous number of lines is the 
method of “k-distribution tabl es’0 (Goody & Yund 119891 ; 
iLacis & OinasllT9^ lEu & Lioulll^2 ). The essence of the 
meth od is to perform Leb esgue, instead of Riemann, integra¬ 
tion (lPierrehumberlll2M^ . when integrating over the opacity 
function of the atmosphere to determine if it is transparent or 
opaque within a given spectral window. Instead of integrating 
over the opacity function itself, which is computationally un¬ 
wieldy as it is hardly a smooth and predictable function, one 
recasts it into its cumulative counterpart—a smooth, mono- 
tonically increasing and computationally pleasing function. 
This cumulative function may then be used to compute the 
transmission function; it is the fraction of radiation passing 
from one layer of the atmosphere to the next within a given 
spectral window. Pigure[T] shows an example of this process. 

The cumulative counterpart of the opacity function is 
known as the “k-distribution function”. The term “k- 
distribution table” is commonly used, because this cumula- 

^ We regard this term as being a synonym, since we will always denote 
opacities hy k and not “k”, following the convention in some paifs of the 
astrophysics literature. However, to preserve tradition we will retain the name 
“k-distribution”. 
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tive function may be tabulated beforehand and then used to 
perfo rm integrations in forward models of radiative transfer 
(e.g. iMarlev et al.lll996t iBurrows et al.lll997t iFortnev et all 
I 2 OIOI) . retrieval models le.g.. iLee. Fletcher & IrwiiT 2012^ 
and t hree-dimensional simu lations of atmospheric circulation 
(e.g.. lShowman et ^l2009^ . 

Nevertheless, several physical and computational issues re¬ 
main either unelucidated or poorly elucidated within the liter¬ 
ature, which provide the motivation behind the current study. 
Our main, physical conclusion is that physical (and not com¬ 
putational) uncertainties associated with the wings of spectral 
lines dominate the error budget. Our technical contribution is 
an ultrafast, open-source computer code to compute the opac¬ 
ity function using modern computing methods and architec¬ 
tures. 

2. METHOD 


T=1500 K, atm, W„=1000 




Fig. 1.— To highlight the salient features explored in this study, we di¬ 
vide the opacity function of the water molecule, as provided by the HITEMP 
database, into three different regions, which we term “bin 1”, “bin 2” and 
“bin 3” in this montage of figures. These bins cover the infrared, infrared- 
optical transitional and optical range of wavelengths. Dividing the opacity 
function into more bins does not alter our qualitative conclusions. As an il¬ 
lustration, we adopt numbers representative of hot exoplanetary atmospheres: 
T = 1500 K and P = 1 atm = 0.98692 bar. Top panel: opacity function 
using spectroscopic quantities from the HITEMP database. Shown are cal¬ 
culations using the full Voigt function and with an ad hoc line-wing cutoff 
of SOOFl. Middle panel: k-distribution functions for the three wavenumber 
different regions of the opacity function. Bottom panel: ttansmission func¬ 
tion corresponding to the three wavenumber regions, both with and without 
the Voigt line-wing cutoff. 


2.1. Theory of k-distribution method versus correlated-k 
approximation 

2.1.1. Restatement of basic theory of k-distributions 

Consider an arbitrary function /(x), where x is the 
wavenumbeiQ normalized by the entire range considered. We 
wish to evaluate the integral over the range Xmin < x < Xmax, 

I = I fix) dx. (1) 

'X Slinin 

Imagine that /(x) may be recast as f{y) such that the quantity 
y is the fractional area under the curve that satisfies / (x) < 
/o, where /o is an arbitrary value of the function. Then, the 
same integral may be evaluated as 

I=f"fiy)dy. (2) 

Jo 

Practically all of the functions we encounter in astrophysics 
may be integrated using this alternative expression. 

More generally, we have 

Fdx = j HdF, (3) 

where TL is the fractional cumulative distribution function of 
another arbitrary function, F(/(x)), that satisfies F < Fq 
and Fq is an arbitrary value of F. In other words, TL gives the 
fractional area under the curve corresponding to F < Fq. 

We make these concepts less abstract by applying them to 
an atmosphere. In the simplest case, suppose that F = f = k, 
where k is the opacity function (with units of cm^ g~^)- It 
follows that f TL dF = J TLk dx = J k dy. The quantity 
y = Jq TL dx is the cumulative sum of intervals. As expected, 
one gets the same answer whether one evaluates f k dx or 
f K dy. A more useful example considers 

f = K, F = exp [—Km), (4) 

where m is the column mass, since the transmission function, 

pOO 

T= / Fdx= F dy^ (5) 

^0 Jo 

^ We use wavenumber, instead of wavelength, because it is the preferred 
choice of spectroscopic databases like HITRAN and HITEMP and spectral 
lines are more evenly spaced across wavenumber (or frequency) than wave¬ 
length. 
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is a quantity that is indispensible for com puting synthetic 
spectra (e.g.. iHeng, Mendonca & Le^ 120141) . The transmis¬ 
sion function is commonly integrated over some wavelength 
range and is the degree or transparency (or opaqueness) of 
this spectral window. For example, in a purely absorbing at¬ 
mosphere the flux passing from one layer to another is given 

by Flayer — F’previous^ “f TtF? (1 7^)? where F’previous tS 

the fl ux from the previous layer and B is the Planck function 
(e.g.. |Heng. Mendonca & Leell2014ll . The second equality in 
equation obtains from expressing the cumulative sum of 
intervals as 

r 

y= 'Hfh dn. (6) 

Jo 

We will refer to K{y) as the “k-distribution function”. 


2.1.2. Correlated-k approximation 

The k-distribution method is exact for a homogeneous at¬ 
mosphere, which almost never happens in practice. For an in¬ 
homogeneous atmosphere, the opacity changes with the tem¬ 
perature and pressure and we have 


r 



K (x) dm 


dx 


K (y) dm 


dy. 


(7) 


making the “correlated-k approximation” and the transmis¬ 
sion function may then be computed as a single integral across 
y. It is the assumption that each value of the cumulative opac¬ 
ity function is always drawn from the same wavenumber in¬ 
terval. 

The mathematics behind the reasoning is identical in the 
case of applying the correlated-k approximation to a homo¬ 
geneous atmosphere with multiple atoms or molecules. For 
illustration, consider only two molecules and a single value 
of the column mass. Let the mixing ratios (relative abundance 
by number) of the molecules be Xi and X 2 . We then have 

T= exp [-X 1 K 1 (yi) rh - ^ 2/^2 ivi) m] dyi 
Jo 

+ / exp [-Xi/ci (y2)m - ^2/^2 (2/2)^] dy2 (H) 
^0 

^2 / exp [-X 1 K 1 (y) rh - X 2 K 2 (y) ni] dy. 

Jo 

Here, the fact that we have two integrals comes from having 
F = exp [—{XiKi + X 2 K 2 )m\ and 

dF =—Fm{XidKi + X2dK2) ■ (12) 

Also, we have 

dyi — 'HfhXi dni ,dy 2 = 'HrhX 2 dK 2 - (13) 


That the k-distribution method cannot be used for an inho¬ 
mogeneous atmosphere may be illustrated using the example 
of a two-layered atmosphere. Each layer has its own opacity 
function and column mass (subscripted by “1” and “2”) and 
the transmission function is 


T= exp [-Ki (yi) mi - K 2 (yi) ^ 2 ] dyi 
Jo 

+ / exp [-Ki (y 2 ) mi - K 2 (y 2 ) ^^ 2 ] dy 2 (8) 
^0 

^2 / exp [-Ki (y) mi - K 2 (y) m 2 ] dy. 

Jo 

That there are two integrals originates from having F = 
exp {—Kiifii — K 2 m 2 ) and 


dF = —F {rhidKi + m2dK2) ■ 


(9) 


Also, we have 


dyi = Hmi dKi ,dy 2 = 'Hfri2 dti 2 - (10) 

The non-equality in equation ® derives from the fact that 
even identical ranges of values in yi and y 2 generally cor¬ 
respond to different ranges of wavenumbers. For example, 
Ki(t/i) and ^ 2 ( 2 / 2 ) are cumulative functions constructed from 
their own cumulative sum of intervals. By contrast, ^ 1 ( 2 / 2 ) 
and ^ 2 ( 2 / 1 ) are cumulative functions constructed from the cu¬ 
mulative sum of intervals of their counterparts, meaning that 
the contributions are drawn from different wavenumber inter¬ 
vals even at the same value of the cumulative sum of intervals. 
Generally, we expect these four cumulative functions to have 
different functional forms. This peculiar property is an un¬ 
avoidable consequence of working with cumulative functions. 

Physically, in employing the k-distribution method, the 
price being paid is that the wavenumber information has been 
scrambled. If one assumes that y = yi = j/ 2 , then one is 


We have intentionally written things out explicitly to illus¬ 
trate the fact that one can avoid dealing with two integrals 
if a single, total opacity function is constructed first (k = 
XiKi -f X 2 K 2 ) before its cumulative function is computed. 

Again, unless yi = y 2 , the two integrals cannot be com¬ 
bined. Since this reasoning holds for multiple molecules 
in a homogeneous atmosphere, it must also hold for multi¬ 
ple molecules in an inhomogeneous atmosphere. We con¬ 
clude that one needs to first add the opacities of the vari¬ 
ous molecules in an atmosphere, weighted by their relative 
abundances, prior to constructing the cumulative function of 
the opacity. If one adds the cumulative opacity functions 
of different molecules, then one is effectively employing the 
correlated-k approximation. 

Both lines of reasoning can be straightforwardly gener¬ 
alised to an inhomogeneous atmosphere containing a single 
atom or molecule and with N layers, a homogenous atmo¬ 
sphere with N atomic or molecular species, or an inhomo¬ 
geneous atmosphere with an arbitrary number of layers and 
species. 

A common source of confusion in the literature is the fail¬ 
ure to distinguish the method (k-distribution) from the ap¬ 
proximation (correlated-k). For example, the “correlated-k 
method” is a misnomer. 

2.2. Implementing the k-distribution method 

Consider equal intervals in x and let the interval be denoted 
by 5x. Such a uniform grid in x generally leads to a non- 
uniform grid in k { x ). Its virtue is that it reduces our problem 
to one of sorting and ordering, since every value of k ( x ) is as¬ 
sociated with 5x (and we do not have to keep track of chang¬ 
ing values of the interval). For a fixed value of the opacity 
(ko)j we count the number of points that satisfy k { x ) < kq . 
If Nx points are counted, then we have 

Wx 6x 


( 14 ) 
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where Aa; = Xmax — a^min is the range of x being considered. 
We also have 5x = Ax/N^, where is the total number of 
intervals in x. It implies that the interval in y is also equal. 


Sy = 


6x 

Ax 


1 


(15) 


By running through all possible values of kq, one constructs 
K(y). Since k(j/) is a monotonic function that is typically 
smoother than k{x), it may be resampled and defined over a 
much smaller number of points, Ny -C N,^. It is then used to 
calculate T for any value of to. 


2.3. Using the HITRAN and HITEMP databases 

The opacity function is a product of two quantities: the in- 
tegrated line strength (S) and the line profile or shape (4)) 
(IGoodv & Yung|[l9^ . 

K = S<^. (16) 


The integrated line strength depends only on the temperature 
(T), while $ depends on both temperature and pressure (P). 
Note that some references collectively refer to opacities (with 
units of cm^ cross sections (with units of cm^) and ab¬ 
sorption coefficients (with units of cm~^) as “absorpt ion co¬ 
efficients” (e.g.. Appendix 2 of iGoodv & Yung|ll989h . Only 
when K is an actual opacity is S' = S(T) with no dependence 
on pressure. 

By invoking the principle of detailed ba lance and lo¬ 
cal thermodynamic equilibrium, one obtains (iPennerl [19521 
iRothman et al.lll99^ . 


92^21 

Sirciy^mQ 



AE\ 

k^J 


1 — exp 


hcu 


(17) 

where p 2 is the statistical weight of the upper level (of a given 
line transition), A 21 is the Einstein A-coefficient, c is the 
speed of light, 1 / is the wavenumber, to is the mean molecular 
mass, Q is the partition function, AE is the energy difference 
associated with the line transition, is Boltzmann’s constant 
and h is Planck’s constant. The partition function relates the 
number density associated with an energy level with the total 
number density and is a function of T. 

In practi ce, a more useful expr ession for the integrated line 
strength is (iRothman et al.|[T 99 CT . 


S Qo / AE AE A 1 “ exp {—hcv/k'QT) 

So Q V kBTo ) I - exp {-hcv/kBTo) ’ 

(18) 

where all of the quantities subscripted with a “ 0” are specified 
at a re ference temper ature, Tp. The HITRA N (IRothman et al.l 
I 2 OI 3 I) and HITEMP (IRothman et al.ll2M^ databases provide 
tabulated values of all of the quantities needed to construct S 
using To = 296 K. 

The Voigt profile is the convoluti on of the Lorentz and the 
Doppler profiles (e.g., iDrainellWIh . 


$ = 




i7v = - 

TT 



exp 

[u — u')^ + a? 


du', 


(19) 


where Td = r'o'\/21n2 k^T is the half-width at 

half-maximum of the Doppler profile, vp is the line-center 


wavenumber, a = s/Io^Tb/Ty) is the damping parameter and 
u = s/\n2{v—vp)/Tjp- Our definitions for Ed, a and M depart 
slightly from the tradit i onal o nes in order to be consistent with 
iLetch worth & Benner|(l2007l) . We have included the effects of 
pressure broadenin g within our definition of the half-wi dth of 
the Lorentz profile lMihalasl[r970HRothman et al.lll996l) . 


rL = 


A21 

/'T\ 

Ticoll 

Q^air {P 

^self) 1 O^self^self 

dnc 

KTp) 


Po 

Po 


( 20 ) 

where the first term after the equality is typically subdom- 
inant. Pressure broad ening is included via an empirical fit 
(IRothman et al.l[T99^ . whose fitting parameters (ncoii, ctair 
and aseif) given by HITRAN and HITEMP. The refer¬ 
ence pressure is Pq = 1 atm = 0.98692 bar. The subscripts 
“air” and “self’ represent air- and self-broadening, respec¬ 
tively. Eor illustration, we assume that they are present in 
equal proportions (Tseif = 0.5P). We also account for a 
pressure-induced shift ((5shift) of the central wavenumber. 


vp^ yp + 


*^shift.P 

Pp 


( 21 ) 


where ^shift is again a tabulated quantity in HITRAN and 
HITEMP. The data for Jshift is usually sparse. 


2.4. Computing the Voigt profile and the line-wing cutoff 
problem 

There are two challenges associated with the Voigt profile. 
The first challenge is computational: it is difficult to eval¬ 
uate efficiently as it is an indefinite integral. Eurthermore, 
we have to compute the Voigt profile multiple times for ev¬ 
ery line and there is an enormous number of lines. To this 
end, we impleme nt Algorithm 916, whic h was originally writ¬ 
ten for MATLAB (IZaghloul & Alill2012h . The essence of the 
algorithm is to first recast Hy as the real part of the (com¬ 
plex) Eaddeeva function and proceed to express it in terms of 
cosines, sines, a scaled complementary error function and sev¬ 
eral series expansions, as stated in equations (13), (15), (16) 
and (17) of iZaghloul & Alii ([TOI^)- The exponential terms in 
the series ex pansions are the bottle neck in terms of compu¬ 
tational cost: IZaghloul & All ([2012) optimise this process by 
combining the three series evaluations within a single loop. 

It turns out that Algorithm 916 is efficient only for small 
values of a and u. Eor a? -\- > 100, we implement 

third-order Gauss-Hermi te quadrature to compute H y as 
stated in equation ( 8 ) of iLetchworth & Benneii (120071) . For 

-f > 106 , we switch from third- to first-o rder Gauss- 
Hermite quadrature (iLetchvvorth & Benneill^OOTl) . Table 1 of 
ILetchworth & Benneil (120071) provides more details on the in¬ 
tegration methods used as a function of a-\u\ space. Our cri¬ 
teria for switching between the three computatio nal methods 
is loosely based on ILetchworth & Benneii OOOTl) and verified 
by testing and trial-and-error. 

The second challenge is physical: the Lorentz, and 
hence the Voigt, profile over-estimates the far wings 
of the line profile due to pre ssure broadening (see 
iFreedman, Marlev & LoddersI 120081 and references therein). 
Even what “far” actually means is not well understood. 
Although this issue dominates the error budget, it is ei¬ 
ther treated as a n ad hoc cutoff ( i n wavenumber) in the 
line wings (e.g., ISharp & BurrowsI 120071: lAmundsen et alJ 
I 2 OI 4 I) . described qualitatively as a problem with no explicit 
cutoff being specified (e.g., IFreedman, Marlev & Lodd^ 
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Fig. 2.— Elucidating the effects of resampling and spectral resolution (Nu). L 
top, middle and bottom rows are for P = 0.01, 0.1 and 1 atm, respectively. For 
due to overfitting. 



m [g cm 

column: k-distribution functions. Right column: transmission function. The 
y = 100, resampling with 20 Chebyshev coefficients results in discrepancies 


20081) or simply left un mentioned (e .g.. lirwin et alJ 120081: 
Madhusudhan & Seageil 120091: iLee. Fletcher & IrwinI 


201 2t B enneke & Seageil 2012t iBarstow et alJ 12013t 
Lee, Heng & IrwinI I2013t [Line et alJ I2013h . It is our 


hope that this issue will be acknowledged more explicitly 
and transparently in future studies involving atmospheric 
radiative transfer and retrieval. 

In the current study, we do not attempt to solve this physics 
problem, which requires a detailed quantum mechanical cal¬ 


culation. In the absence of a complete, first-principles the¬ 
ory, we instead compare calculations with the full Voigt pro¬ 
file versus those with some arbitrary line-wing cutoff spec¬ 
ified, which we nominally take to be 500 Lorentz widths. 
We emphasize that there is no sound physical reason behind 
choosing this particular cutoff. It is merely used as a proof- 
of-concept comparison against calculations utilizing the full 
Voigt profile. 
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2.5. GPU computing: memory types and parallelization 

We develop our custom-built code (HELIOS-K) using the 
native language of the NVIDIA GPUs, CUDA (Compute Uni¬ 
fied Device Architecture), which is basic ally an embellished 
versio n of the C programming language (I Sanders & KandroB 
[20I3). a major advantage provided by a GPU is the large 
number of computational cores per card 1000 ) for a very 
low cost $1 per core). When compared head-to-head, a 
single GPU will always lose out against a single CPU in terms 
of both computational power and memory—the point is that 
one wins by throwing many, many more GPU cores at the 
problem. A set of 32 consecutive threads is called a “warp” 
and it is crucial that every warp performs exactly the same op¬ 
eration in order to optimise performance. If not, a “branch di¬ 
vergence” occurs and some operations are performed in serial 
operation mode. Each calculation is performed on a thread 
and all of the threads are organised into blocks. 

An indispensible part of writing ultrafast CUDA code is to 
understand the memory design and types on a GPU. Global or 
device memory is the most abundant and can be accessed by 
every thread, but is generally the slowest type. Shared mem¬ 
ory is faster, but may only be accessed by threads within the 
same block. Typically, a well-written CUDA kernel (usually 
called a “function” in other languages) reads data from global 
into shared memory, performs the necessary arithmetic oper¬ 
ations and writes back to global memory. Another bottleneck 
is the passing of information (communication) between the 
CPU and GPU. Exploiting the order-of-magnitude speed-ups 
a GPU has to offer is an exercise in shrewd memory and com¬ 
munication management. Rather than describe each and every 
computing trick we used, we highlight the main ones and refer 
the reader to our open-source code. 

Eor our application, we need to compute the Voigt profile 
for an enormous number of spectral lines across an even larger 
number of grid points in wavenumber. Eurthermore, we need 
to repeat this calculation for multiple combinations of tem¬ 
perature and pressure. It is impossible to perform this com¬ 
putation in a single step, but we may perform a serial loop 
across the lines and parallelise across wavenumber. This al¬ 
lows us to accumulate values of k{x) directly within a register 
(i.e., fastest available memory) without additional write-outs 
to global memory. 

2.6. Sorting and resampling 

Parallel sorting on a GPU is a non-trivial task. Eortunately, 
this has already been implemented as part of the CUDA library 
(https : / /developer . nvidia . com/Thrust). Once 
we have computed k{x), the challenge is to perform the sort¬ 
ing within each bin. Each bin has a width Aa; and the number 
of bins typically used is ^ 10-10^. Sorting each bin in a 
serial fashion would be inefficient when the number of bins 
becomes large. Instead, we sort the entire opacity function all 
at once, but keep track of the bin number each opacity point 
belongs to, which ultimately allows us to reconstruct k(j/) in 
the individual bins. 

Once we have sorted k{x) and obtained K{y), we wish to 
resample K{y) such that it is defined using a much smaller 
number of points (by orders of magnitude). Numerous re¬ 
sampling strategies exist, including least-squares fitting, fast 
Eourier transforms, etc. We find that using a least-squares fit 
with Chebyshev polynomials gives the best outcome in terms 
of accuracy and efficiency, especially since one may exploit 
the recurrence relations to generate Chebyshev polynomials 


of different orders. We perform the fit on In K{y) to avoid nu¬ 
merical oscillations. The least-squares fitting essentially in¬ 
volves solving AC = D for the vector of Chebyshev coeffi¬ 
cients (C), where D is the data vector. Directly computing the 
inverse of the matrix A is expensive; instead, we implement 
“Q-R decomposition” to obtain C (iPress et alJl2007l) . The fi¬ 
nal product of this step is a set of 20 Chebyshev coefficients 
describing K{y) for each bin. 

3. RESULTS 

Unless otherwise stated, our results are based on comput¬ 
ing a pure-water opacity function using the HITEMP line list, 
which consists of ^ 10® spectral lines of water. We empha¬ 
size that this is a proof of concept and that HELIOS-K may 
be used for general mixtures of atoms and molecules. 

3.1. Basic setup 

We base the discussion of our results on a fiducial setup. 
We focus on computing the opacity function for the wa¬ 
ter molecule, since it has the most lines among the major 
molecules expected (compared to CO and CO 2 ) and has the 
least controversial line list available (compared to CH 4 ). In 
Pigure[T] we show two instances of the opacity function; one 
computed using the full Voigt profile and the other with a 
line-wing cutoff applied. We divide the wavenumber region 
into three equal ranges: 0.5-8573.5 cm“^ (infrared to near- 
infrared; > 1.2 ^m), 8573.5-17146.5 cm“^ (near-infrared 
to optical; 0.6-1.2 p,m) and 17146.5-25719.5 cm“^ (optical; 
0.4-0 .6 ^m). Each bin has a width of ISu = 8573 cm“^. 
Within each bin, we adopt a resolution of = ^vjSv = 
10 ®; we will demonstrate later that this attains convergence. 
Our results point to the same qualitative conclusions even 
when more bins are used (not shown). 

It is readily apparent that the choice of cutoff is a significant 
source of error in the near-infrared and optical, because it af¬ 
fects the weak lines more strongly, even prior to the mapping 
of the opacity function to its k-distribution counterpart. We 
emphasize that this problem remains, regardless of whether 
the k-distribution method is used. Eor the k-distribution func¬ 
tion and the transmission function, the influence of the choice 
of line-wing cutoff is seen to be significant. We will investi¬ 
gate this issue in more detail. 

3.2. Resampling as an insignificant source of error 

A necessary, intermediate step to check is whether the re¬ 
sampling of the k-distribution function using least-squares fit¬ 
ting introduces a significant source of error to our results. In 
Eigure|2] we compute the transmission function in two ways; 
using the direct output from the mapping of k{x) to K{y) and 
the resampled K{y). The difference between the two calcu¬ 
lations is typically <C 1% when N^, > 10®. Remarkably, 
resampling is not a significant source of error independent of 
the value of the column mass, i.e., it is equally robust in both 
optically thin and thick parts of the atmosphere. 

3.3. Choosing the correct bin resolution 

Even though convergence within each bin is tied to the 
number of lines present, we find that an easier rule of thumb 
is to use a minimum value of as a convergence criterion. 
Eigure |2] shows that convergence is comfortably attained for 
> 1000. This conclusion holds even when 1000 bins are 
used (not shown). 
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m [g cm“^] 

Fig. 3.— Transmission function subjected to different choices of the line¬ 
wing cutoff and at different pressures. The cuts of 100 cm~^ and 500 Fl are 
chosen to match each other at P = 1 atm. 


T=1500 K, P=1 atm, W =1000, 1000 bins 



Fig. 4.— Fractional difference in transmission function, between calcula¬ 
tions using the full Voigt profiles and with a cutoff of Fl imposed, for 1000 
bins across the same wavenumber range. 

3.4. Line-wing cutojfas the largest source of error 

We further explore our claim that the line-wing cutoff is 
the largest source of error in computing and using an opac¬ 
ity function, regardless of whether one uses the k-distribution 
method. In Figure [3l we show different calculations of T for 
various cutoff choices: an absolute cutoff (of 100 cm~^, fol¬ 
lowing the choice made bv ISharp & Burrowsl2007^ and an ad 
hoc cutoff of 500 Lorentz widths. These choices are made 
such that they produce the same results at P = 1 atm. At 
higher pressures, we see that deviations appear. For a given 
value of the column mass, the error is ~ 10% to even a factor 
of several in some instances. 

We quantify this error in more detail. Figure |4] shows the 
fractional difference in T, between calculations using the full 
Voigt profiles and those with a cutoff of Fl imposed, for 1000 
bins across the same wavenumber range. Across a broad 
range of column masses (10“^ < m < 10^ g cm“^), we 
compute the median, mean and maximum fractional differ¬ 
ences using the full-Voigt calculations as a baseline compar¬ 
ison. (We emphasize this does not imply that using the full 
Voigt profile is correct.) The median and maximum fractional 
differences are dominated by small and large column masses, 
respectively, and are not representative, but we show them 



line wing cutoff [cm 

Fig. 5.— Performance of HELIOS-K broken down into the various compu¬ 
tational tasks: reading in line-list data (“Input”), computing Lorentz width, 
Doppler width, line strength and line center shift (“Line”), computing Voigt 
profiles and the opacity function (“«;(i/)”), sorting values of k{x) into K{y) 
(“Sort”), resampling K{y) (“Resampling”) and computing T for 1000 val¬ 
ues of rh (“Transmission”). For this plot only, we are using the HITRAN 
database with ~ 10® water lines. 

for completeness. We see that the mean fractional difference 
is ~ 10% across all wavenumbers for P = 1 atm, imply¬ 
ing that a similar uncertainty is present for the computed flux 
or synthetic spectrum. We expect that the median fractional 
differences are larger for higher pressures. Elucidating the 
full consequences of the uncertainty, associated with pressure 
broadening, for calculations of radiative transfer and retrieval 
is deferred to future work. 

Generally, we And that the uncertainties associated with 
the line-wing cutoff are typically larger than those due to, 
e.g., resampling, as long as a sufficient bin resolution is used 
{Ny > 10^, as previously demonstrated). 

3.5. Performance 

We execute performance tests on a NVIDIA Tesla K20 
GPU card, which has 2496 cores. For these tests, we use 
the FIITRAN 10^ water lines), instead of the HITEMP, 
line list, as the entire calculation fits within a single K20 
GPU card. (The HITEMP water line list is provided in 34 
separate chunks, which we simply load in serial.) Eigure |5] 
breaks down the performance of our code, which we name 
HELIOS-K, in terms of the various tasks executed. Unsur¬ 
prisingly, the computational cost goes up with bin resolution 
and line-wing cutoff. Generally, HELIOS-K takes ^ 1 s to 
compute ^ 10^ spectral lines of water. We anticipate that such 
a level of performance allows for efficient and broad sweeps 
of the parameter space of exoplanetary atmospheres. 

4. DISCUSSION & CONCLUSIONS 

4.1. Towards uniform standards: a checklist for opacity 
function calculations 

The details of how opacity functions are computed and used 
by various studies in the literature remain vague or incom¬ 
plete. We suggest that a path towards uniform standards in¬ 
volves explicitly addressing the following questions (and pub¬ 
lishing the answers to them). 

• Does the study claim a “line-by-line” calculation of the 
opacity function (e.g., iMadhusudhan & Seaged 1^091: 
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iBenneke & Seage3l2012l) ? If so, are the lines being 
sampled in an adequate way? E.g., if there are Mines 
lines, is Mampie > Mines, where Mampie IS the num¬ 
ber of wavenumber/wavelength points used? If spe¬ 
cial circumstances (e.g., very broad lines) allow for 
Mampie ~ Mines to be justified, has this been demon¬ 
strated explicitly? Does the study show results from 
convergence tests? Often, what are effectively opacity¬ 
sampling techniques (Mampie ^ Mines) are mislead¬ 
ingly claimed as being “line-by-line”. 

• How is the Voigt profile being computed? Is it being 
directly evaluated as an indefinite integral? Or has a 
transformation and/or approximation(s) been taken? 

• If the k-distribution method is adopted, how many bins 
are specified? How is the opacity function resampled 
within each bin, i.e., what is the resampling method? 
Has the study demonstrated that an adequate intra-bin 
resolution has been used? 

• Are k-distribution tables separately computed for each 
molecular species and then added together—weighted 
by the relative abundance of each species—afterwards? 
If so, then the correlated-k approximation has been used 
and this should be explicitly mentioned. 

• Is pressure broadening being considered? If so, is the 
study imposing line-wing cutoffs? Is the cutoff speci¬ 
fied as an absolute number or as a specific number of 
Lorentz or Doppler widths? Have the uncertainties as¬ 
sociated with this choice been explored and quantified? 
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The preceding checklist may be a useful guide for review¬ 
ing studies that perform radiative transfer or retrieval calcula¬ 
tions. 


4.2. Summary 

We have constructed an open-source, ultrafast, GPU code 
written using CUDA, named HELIOS-K, which takes a line 
list as an input and computes the opacity function of the at¬ 
mosphere for any mixture of atoms and molecules. The dom¬ 
inant source of error stems from an unsolved physics prob¬ 
lem: describing the far line wings of spectral lines affected 
by pressure broadening. In the absence of a complete theory, 
we (and others before us) have applied an ad hoc cutoff of 
the line wing for our calculations. Notwithstanding this is¬ 
sue, HELIOS-K provides the exoplanet community with an 
efficient tool for computing opacity functions. 
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